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Abstract 

Based on an analysis of the short range chemical environment of each atom in a system, standard 
machine learning based approaches to the construction of interatomic potentials aim at determining 
directly the central quantity which is the total energy. This prevents for instance an accurate 
description of the energetics of systems where long range charge transfer is important as well as 
of ionized systems. We propose therefore not to target directly with machine learning methods 
the total energy but an intermediate physical quantity namely the charge density, which then in 
turn allows to determine the total energy. By allowing the electronic charge to distribute itself 
in an optimal way over the system, we can describe not only neutral but also ionized systems 
with unprecedented accuracy. We demonstrate the power of our approach for both neutral and 
ionized NaCl clusters where charge redistribution plays a decisive role for the energetics. We are 
able to obtain chemical accuracy, i.e. errors of less than a milli Hartree per atom compared to the 
reference density functional results. The introduction of physically motivated quantities which are 
determined by the short range atomic environment via a neural network leads also to an increased 
stability of the machine learning process and transferability of the potential. 
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Atomistic simulations for materials are nowadays widely applied to understand and de¬ 
sign materials. A wide range of simulation methods exist, ranging from quasi exact many 
electron wavefunction methods [I], over density functional theory (DFT) calculations [2] and 
semi-empirical quantum mechanical methods to inter-atomic potentials [5j. For the 

different methods there are well known trade-offs between the computational costs and the 
accuracy of the calculations. In contrast to classical force fields, density functional calcula¬ 
tions give sufficient accuracy for a wide range of properties and are therefore the method of 
choice in a huge number of studies. Due to their computational cost, DFT simulations are 
however limited in practice to systems containing less than a few thousand atoms. Since 
there is a great need to do highly accurate simulations for larger system there have been 
numerous efforts to improve the accuracy of force fields, without increasing their cost too 
much. It has been widely recognized that fixed charges limit the accuracy of the estab¬ 
lished standard force fields. In the chemistry and biology community polarizable force fields 
have therefore been developed [5]. However such polarizable force fields allow only for the 
displacement of charges, but do not allow for a true charge transfer over large distances. 
Charge equilibration methods m offer this possibility and have been implemented in force 
fields such as ReaxFF j9]. 

For covalent bulk materials such as carbon and silicon machine learning based total 
energy schemes mm turned out to give density functional accuracy at greatly reduced 
numerical cost. For finite systems such as clusters, the construction of highly accurate 
machine learning based force fields is more difficult, since atoms at a surface behave quite 
differently from atoms in the bulk. It is not surprising that an analysis of the charge 
distribution of NaCl clusters obtained from density functional calculations clearly reveals 
a charge transfer between atoms at the surface and in the core of the cluster. The fixed 
charges of ±1 electron used in the established force fields such as the Tosi-Fumi force field [5] 
are therefore clearly inadequate to describe such systems with high accuracy. Since in the 
standard machine learning based interatomic potentials the energy of the whole system is 
written as a sum of the energies of individual atoms na. and since the energy of an individual 
atom is determined by the short range chemical environment, long range charge transfers as 
well as ionization can not well be described. 

Constant charges are also present in the majority of force fields for bio-molecules and 
finding good charges for the use in such force fields is highly non-trivial. Even though the 
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electrostatic part of the interaction energy is not the dominating one in these force fields, it is 
to be expected that variable charges could also lead in this context to considerable improve¬ 
ments of the accuracy. Environment dependent charges obtained from neural networks (NN) 
were recently introduced Pi, but this approach does not give the correct total charge of the 
system and can for instance lead to fluctuations of the total charge in molecular dynamics 
simulations. The fact that the total charge can not be fixed also prevents the treatment of 
ionized systems. In addition, atomic reference charges must be provided in this approach. 
The extraction of such charges from ab initio calculations is always ambiguous to a certain 
extent. 

A machine learning based approach, which in the spirit of our paper, does not directly aim 
at the total energy, but at an intermediate physical quantity is the work of Snyder IHII5] 
and coworkers, where they propose to construct machine learning based kinetic energy func¬ 
tionals. Whereas in that work the charge density can in principle adopt any form we restrict 
our charge density to be an approximate superposition of atomic charge densities. In this 
way we exploit the well known fact that the charge density in molecules, clusters and solids 
is given within a first approximation by a superposition of atomic charge densities. The dis¬ 
tribution of the electronic charge density is determined by atomic environment dependent 
electronegativities. These electronegativities are in turn determined by the short range envi¬ 
ronment of the associated atom and can easily be predicted by a neural network process. To 
demonstrate the power of our approach we choose ionic clusters where a correct description 
of the charge density is essential since bonding is mediated through charge transfer. The 
fact that the charge distribution can redistribute itself over the whole cluster will allow us 
to treat both surface and bulk atoms with very high accuracy and we will demonstrate that 
density functional accuracy can be obtained with such an approach for clusters. In contrast 
to conventional force fields, our approach also allows to describe ionized systems without 
the need of any reparametrization. 

We consider a system consisting of y Sodium (Na) and y Chloride (Cl) atoms. We 
postulate the following form for the total energy expression 
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FIG. 1. Schematic illustration of the proposed method. {R4 and {Gj} are atomic positions and 
symmetry functions, respectively. 

where E? t are some fixed reference energies which in our implementation are set to energies of 
isolated atoms, qi's are atomic charges, Xt is the environment dependent atomic electroneg¬ 
ativity of atom i whose functional dependence is determined by a neural network approach. 
To describe the hardness [16jj of atom i we introduce element dependent atomic hardness 
terms The charge density of the system, p(r), is in our approach assumed to be a su¬ 
perposition of spherically symmetry Gaussian functions centered at atomic positions, each 
normalized to the corresponding atomic charge qi given by 
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With this choice for the atomic charge densities, the total energy of Eq. (|Tj) can be calculated 
analytically. 
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where 7 i7 - = , l n and r„- the distance between atoms i and j. The atomic charges q % are 
V a i+ a j 

implicitly environment dependent through Xi as w iH be explained below and are therefore 
implicit functions of the atomic positions. The implicit dependence of the q^s is obtained 
by minimizing the total energy of Eq. ([ 2 ]) with respect to the g,’s. This leads to a linear 
system of equations where all matrix elements depend on the atomic positions, 
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The fact that the electrostatic interaction energy of any continuous charge density is always 
positive and that the J^s are positive constants, implies that the matrix A is positive definite 
and that the linear system of equations has therefore always a unique solution which gives 
the minimal electrostatic energy. The linear system Eq. ([3]) has to be solved under the 
constraint that the atomic charges sum up to the correct overall charge q to t = YliLiQi- 
Adding this constraint via the Lagrange multipliers leads to the modified linear system of 
equations 


AQ = X, (4) 

where Q and X are vectors of the dimension (N + 1) and A is a (N + 1) x (N + 1) matrix. 
In expanded form Eq. Q is given as, 
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In this way, we can allow for charge transfer over long distances without the need to find 
the presumably extremely complicated explicit long range environment dependence of the 
qi s. All that is needed to get the implicit long range dependence of the q, 's on the atomic 
positions, is an explicit short range dependence of the yds, which is fixed once and for all 
by our neural network together with the solution of a simple linear system of equations. In 
addition, the total charge of the system is conserved unlike the method given by Ref. [13]. 
This approach is physically motivated, since in a Kolin Sham density functional calculation 
the total energy is minimized with respect to the charge density distribution. So the approach 
can be considered as some kind of constrained minimization of the total energy, where the 
form of the charge density is restricted to be a sum over Gaussian functions with amplitudes 
qi. In Kohn Sham density functional theory, the total energy consists of course not only of 
the electrostatic term, but in addition there are the kinetic and exchange correlation terms 
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FIG. 2. (color online) Error distribution of the training and validation sets. In total the training 


and validation sets contain 15191 and 7658 structures, respectively. 



FIG. 3. (color online) Energy versus variation of compression and expansion ratios relative to the 
equilibrium geometry of the respective method for a 64-atom NaCl cubic structure. 

which oppose a charge distribution that would merely minimize the electrostatic part of the 
total energy. In our scheme this opposing force is represented by the constraints on the form 
of the charge density and the environment dependence of the Xi s - 

Fig. [3] shows LDA, PBE [T?J and NN energies for different compression and expansion 
ratios relative to the equilibrium geometry of the respective method for a 64-atom NaCl 
cubic structure. Even though our energy expression is missing the kinetic and exchange 
correlation terms, the physically important energy differences can be calculated in exactly 
the same was as in any density functional scheme based on the Hellmann-Feynman theorem. 
This theorem tells us that energy differences are given by the integral over the atomic forces 
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FIG. 4. (color online) Correlation of PBE/NN energies versus that of LDA for neutral test points. 


times the displacements. The force acting on the atoms is the classical electrostatic force 
arising from the charge distribution obtained through quantum mechanics. Since our charge 
density is determined by a variational principle it obeys the Hellmann-Feynman theorem. 
In this way it is actually also granted that the charges obtained by solving the linear system 
of equations give an accurate charge density within the limitations imposed by our adopted 
form of the charge density. Unreasonable charge densities would lead to wrong forces and 
hence to wrong energy differences between different structures. 

A schematic illustration of our method is given in Fig. [l] Functional forms of symmetry 
functions ({(?*}) are those given in Ref. [TH] with a modification of the cutoff function which 
is in our implementation a polynomial. In order to construct and examine the method, we 
generated a large number of reference data points using the BigDFT [19] code. An accurate 
evaluation of the electrostatic term in density functional calculations is of great importance 
in particular for ionized clusters. In the BigDFT code, the Hartree potential is calculated 
using the method given in Ref. [2D] which enables us to have accurate energetics for the 
ionized reference data points. The DFT calculations are performed with the local density 
approximation (LDA). A large number of data points consist of NaCl neutral and ionized 
(+1 and +2) clusters ranging from 8 to 80 atoms with step sizes of 4 are generated. Data 
points are split into three sets; training, validation and test. Clusters with fewer than 44 
atoms are considered in training set. The structures of 64-atom cluster are reserved for 
testing and the rest are used as validation data set. 
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In order to avoid overfitting, all training runs are performed with only 15 epochs. The 
activation function of the output layer (values of electronegativities) is taken to be the tanh 
function rather than the commonly used linear function. This allows us to avoid a strong 
variation of atomic electronegativities from one atom to another within a structure, which 
in turn may result in too large variations in atomic charges. In contrast to the standard 
NN potentials in Refs. mm, we can achieve in our approach small errors with a small 
number of nodes in NN hidden layers. Among several fits using various NN architectures 
and several training runs with different initial random numbers for NN weights, we obtained 
the best compromise between small RMSE and transferability of the potential for training 
and validation data sets with the architecture 51-3-3-1. Gaussian widths of 1.0 and 2.0 bohr 
are consider for sodium and chlorine atoms, respectively. Values 0.2 and 0.1 in atomic units 
are found to be suitable for .J tl for sodium and chlorine atoms, respectively. 

The obtained RMSE of the training and validation sets is 0.26 mHa per atom, the RMSE 
of the neutral test set is 0.13 mHa per atom and the RMSE of the ionized test set (including 
q tot = +1 and q tot = +2) is 0.44 mHa per atom. The error distribution of the training 
and validation sets is illustrated in Fig. [2] For a few structures, the error is about 1 mHa 
per atom while for the majority of structures, the accuracy of our method is much higher 
than the so-called chemical accuracy of lkcal/mol. Fig. [3] shows LDA, PBE na and NN 
energies for different compression and expansion ratios relative to the equilibrium geometry 
of the respective method for a 64-atom NaCl cubic structure. The result shown in Fig. [3] 
demonstrates clearly the transferability of our method since such compression and expansion 
were neither included in the training nor in the validation set. 

Fig. [4] illustrates the correlation between PBE/NN and LDA energies. All energies in 
Fig. [4] are relative to the average energy per atom of all structures for each method. Corre¬ 
lation between NN and LDA energies is better than that between PBE and LDA both for 
low- and high-energy structures, which shows that in principle even accuracies higher than 
those obtainable from density functional theory could be achieved if the training/validation 
set energies are calculated with a more accurate method. 

In order to check whether the entire low energy configurational space is well described 
by our interatomic potential we performed minima hopping [2T] (MH) runs for different 
cluster sizes. MH explores in a systematic way this part of configurational space and can 
thus detect whether the interatomic potential fails to describe certain regions of this space. 


No physically unreasonable configuration was found in all these runs indicating that the 
interatomic potential describes well the entire low energy configurational space. 

Based on a machine learning algorithm, a novel scheme is presented to reproduce with 
high accuracy potential energy surfaces of quantum mechanical origin. Instead of predict¬ 
ing atomic energies directly, we use NN methods to predict environment dependent atomic 
electronegativities from which atomic charges are obtained by a charge equilibration pro¬ 
cess m- Once the charges are available the total energy can be calculated easily. Applying 
the method to neutral and ionized sodium chloride clusters shows that unprecedented accu¬ 
racy can be obtained for a particularly difficult system, namely clusters where the atomic 
environment differs drastically between surface and core atoms. The error in energies of our 
interatomic potential compared to the reference density functional data is far less than the 
error in total energies arising from the use of different exchange correlation functionals. The 
potential is highly transferable and does not give rise to any unphysical structures. 
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